clear; clc; close all;

%% 1. 参数设置
mu = 0.2;
x0 = 0.2;     % 初始状态
q0 = 0.1;
k_start = 1.6;  % k 的起始
k_end   = 2.1;  % k 的结束
dk      = 0.002;% k 的步长（可调）
k_vals  = k_start : dk : k_end;

% 迭代相关设置
transient_steps = 300;  % 前期过渡步数
record_steps    = 100;  % 用于分岔图记录的后续步数

% 用于存放计算结果
LE1_array = zeros(size(k_vals));
LE2_array = zeros(size(k_vals));

x_bif = [];  % 用于存储分岔图 x
q_bif = [];  % 用于存储分岔图 q
k_bif = [];  % 对应的 k

%% 计算 LEs 和分岔
index = 1;
for k = k_vals
    % 每次从同一个 (x0,q0) 开始
    x = x0;
    q = q0;
    
    % a) 先做 transient_steps 步迭代，排除瞬态
    for i = 1:transient_steps
        [x, q] = MLM(x, q, mu, k);
    end
    
    % b) 开始记录 record_steps 步，用于分岔图
    x_record = zeros(1, record_steps);
    q_record = zeros(1, record_steps);
    x_record(1) = x;
    q_record(1) = q;
    for i = 2:record_steps
        [x_record(i), q_record(i)] = MLM(x_record(i-1), q_record(i-1), mu, k);
    end
    
    % 拼接到分岔图数据中
    x_bif = [x_bif, x_record];
    q_bif = [q_bif, q_record];
    k_bif = [k_bif, k * ones(1, record_steps)];
    
    % c) 计算李指数(LE1, LE2)
    [LE1, LE2] = LEs_MLM(mu, k, x0, q0);
    
    LE1_array(index) = LE1;
    LE2_array(index) = LE2;
    
    index = index + 1;
end

%% 绘图
% 上图：LE1, LE2 - k
subplot(2,1,1);
hold on; grid on;
plot(k_vals, LE1_array, 'r', 'LineWidth',1.5, 'DisplayName','LE_1');
plot(k_vals, LE2_array, 'b', 'LineWidth',1.5, 'DisplayName','LE_2');
yline(0,'k--','LineWidth',1);  % 零线
xlabel('k'); ylabel('LEs');
title(['2D-MLM: \mu=',num2str(mu), ', (x_0,q_0)=(',num2str(x0),', ',num2str(q0),')']);
legend('Location','Best');

% 下图：分岔图 (x vs. k) 与 (q vs. k)
subplot(2,1,2);
hold on; grid on;
% 建议把 x, q 分两条曲线画在同一个坐标系，也可分开
scatter(k_bif, x_bif, 1, 'g', 'filled', 'DisplayName','x');
scatter(k_bif, q_bif, 1, 'm', 'filled', 'DisplayName','q');
xlabel('k'); ylabel('x, q');
legend('Location','Best');

%% ================ 子函数：迭代映射 =================
function [x_next, q_next] = MLM(x, q, mu, k)
    % 2D-Memristive Logistic Map:
    % x_{n+1} = mu*x_n*(1 - x_n) + k*cos(q_n)*x_n
    % q_{n+1} = q_n + x_n
    x_next = mu*x*(1 - x) + k*cos(q)*x;
    q_next = q + x;
end

%% ================ 子函数：计算 2DMLM 的 LEs =================
function [LE1, LE2] = LEs_MLM(mu, k, x0, q0)
    % 使用离散系统的雅可比矩阵 + QR 正交化方法，迭代若干步后求平均指数
    % 可根据需要增加/减少迭代次数
    MaxIter = 1000;
    Transient = 200;  % 舍去前 200 步
    
    x = x0; 
    q = q0;
    % 初始正交基（2D）:
    e1 = [1; 0];
    e2 = [0; 1];
    
    sum_log1 = 0;
    sum_log2 = 0;
    
    for i = 1:MaxIter
        % 1) 计算雅可比矩阵 J(x,q)
        %    x_{n+1} = f1(x,q), q_{n+1} = f2(x,q)
        %    f1 = mu*x*(1-x) + k*cos(q)*x
        %    f2 = q + x
        %
        %    对 x_{n+1} 关于 x,q 求偏导:
        %      df1/dx = mu*(1-2x) + k*cos(q)
        %      df1/dq = -k*sin(q)*x
        %
        %    对 q_{n+1} 关于 x,q 求偏导:
        %      df2/dx = 1
        %      df2/dq = 1
        J = [ mu*(1-2*x) + k*cos(q),   -k*sin(q)*x;
              1,                      1            ];
        
        % 2) 用 J 作用在正交基上
        e1_new = J * e1;
        e2_new = J * e2;
        
        % 3) QR 分解正交化
        [Q, R] = qr([e1_new e2_new],0);  % 2x2
        e1 = Q(:,1);
        e2 = Q(:,2);
        
        % 4) 累加对数范数
        if i > Transient
            sum_log1 = sum_log1 + log(abs(R(1,1)));
            sum_log2 = sum_log2 + log(abs(R(2,2)));
        end
        
        % 5) 做一次迭代更新
        [x, q] = MLM(x, q, mu, k);
    end
    
    % 计算平均指数
    eff_iter = MaxIter - Transient;  % 有效迭代次数
    LE1 = sum_log1 / (eff_iter);
    LE2 = sum_log2 / (eff_iter);
end
